Quasi-non-local gradient-level exchange-correlation approximation for metals and alloys 
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The flexibility of common generalized gradient approximation for the exchange-correlation energy is inves- 
tigated by monitoring the equilibrium volume of transition metals. It is shown that no universal gradient-level 
approximation yielding consistent errors for all metals exists. Based on an element-specific optimization, the 
concept of quasi-non-local gradient-level approximation is introduced. The strength of the scheme is demon- 
strated on several transition metal alloys. 



Density functional theory (DFT)I has become a funda- 
mental first-principles research tool within the modern ma- 
terials science. The great breakthrough during the last 30- 
40 years should primarily be attributed to the local density 
formalism, 2 where the many-body interactions are efficiently 
incorporated within an effective one-electron local potential. 

The first description of the unknown exchange-correlation 
potential was provided by the local density approximation 
(LDA), assuming uniform electron density on the scale of the 
exchange-correlation hole. This seemingly rough approxima- 
tion turned out to be unexpectedly accurate in total energy 
calculations, which in fact made the early DFT success pos- 
sible. Attempts to go beyond LDA within the framework 
of the local density formalism have led to the elaboration 
of the density gradient corrected functionals. Accordingly, 
two major gradient-level density functional (GDF) families 
emerged. The generalized gradient approximation (GGA) 3 ~ 6 
aimed to stabilize the diverging term from the second order 
gradient expansion, 2 and gave for the first time a proper de- 
scription of many important solids, such as the ferromag- 
netic iron. The subsystem functional approach (SFA)£ orig- 
inated from the nearsightedness principle 8 and incorporates 
inhomogeneous electron density effects through well adapted 
model systems. The simplest SFA was put forward within the 
Airy gas approximation, 9 which was later further developed 
into various gradient-level approximations jiSi 2 . For both GDF 
families, LDA represents the lowest order approximation and 
thus the correct limit in systems with densities showing negli- 
gible inhomogeneities. 

A GDF exchange-correlation energy depends on the elec- 
tron density n and its gradient Vn, viz. 

£ x G c DF M = / d 3 rf xc (r s ,s), (1) 

with r s = (3/(4 7 rn)) 1 / 3 and s = \\7n\/ (2n(3TT 2 n) 1 / 3 ). 
The exchange-correlation energy density f xc (r s ,s) is com- 
monly expressed with the help of the enhancement function 



F xc (r s , s) over the exchange energy density of the uniform 
electron gas, e^ DA (n). For slowly varying electron density 
(s -> 0), F xc (r s ,s) F^ A (r s ) = 1 + e^ A (n) / el BA (n) 
where e^ DA (n) is the LDA correlation. The explicit form of 
F xc (r s ,s) contains all information about the actual approx- 
imation. A detailed comparison of different enhancement 
functions from the GGA and SFA families 13 shows that the 
behavior of F xc (r s , s) in the rapidly varying density regime 
(corresponding to s > 0.5 for r s < 4 and s > 0.2 for r s > 4) 
determines the performance of the approximation for a partic- 
ular inhomogeneous electron system. 

In this Letter, first we investigate the flexibility of the 
GDFs by monitoring their behavior for bulk transition met- 
als. Based on the SFA concept, next we introduce a quasi- 
non-local gradient-level approximation (QNA) that performs 
equally well for mono-atomic systems as well as for multi- 
component solid solutions. For sake of transparency, we 
limit the present assessment to two familiar GGA exchange- 
correlation density functionals, namely the Perdew-Burke- 
Ernzerhof (PBE) 4 and the revised PBE (PBEsol) 6 functionals. 
Nevertheless, the main conclusion and the emerging QNA can 
easily be extended to any GDF type of description. 

The PBE and PBEsol approximations are based on sim- 
ilar enhancement functions. The prior was designed to pro- 
vide accurate atomic energies whereas the latter was opti- 
mized for bulk and surface systems by restoring the original 
gradient expansion behavior for the exchange part and adjust- 
ing the correlation term using the jellium surface exchange- 
correlation energies obtained at meta-GGA level. 14 The pa- 
rameters controlling these effects are /i and f3 and their values 
for PBE and PBEsol are {p, /3) PBE = (0.21951,0.066725) 
and 0,/3)pbeso1 = (0.123457, 0.046000), respectively. The 
large (44 and 31%) difference between the two sets of param- 
eters indicates that the atomic and bulk regimes require rather 
different enhancement functions. The LDA functional is re- 
covered for (fi, /3)lda = (0, 0). 

We demonstrate the performance of different GDFs' by 
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TABLE I, Equilibrium Wigner-Seitz radii (in units of Bohr) at K 
for a selected set of transition metals. Experimental values are the 
same as those used in Ref. [23| (Cu, Pd), Ref. 24 (Fe) and Ref. 
l25l (Au, W). For Nb and V, the experimental values were obtained 
by extrapolating the room temperature lattice constants 26 to K us- 
ing the experimental thermal expansion coefficients 27 and the Debye 



temperature i 2 ^ 


Element 


TOPBE 


™PBEsol 


™(expt) 


V 


2.789 


2.752 


2.805 


Fc 


2.641 


2.600 


2.661 


Cu 


2.687 


2.638 


2.661 


Nb 


3.080 


3.040 


3.066 


Pd 


2.923 


2.873 


2.866 


W 


2.970 


2.942 


2.940 


Au 


3.084 


3.028 


3.002 



calculating the equilibrium volumes of seven transition metals 
(V, Fe, Cu, Nb, Pd, W, Au) and four alloys (V-W, V-Fe, CuAu 
and CU3AU). The electronic structure and total energy cal- 
culations were performed using the exact muffin-tin orbitals 
method.—— The Kohn-Sham equations were solved within 
the scalar relativistic approximation, the Green's function was 
calculated for 16 complex energy points distributed exponen- 
tially on a semicircular contour including the valence states 
and employing the double Taylor expansion approach.— The 
basis set included s, p, d, and / orbitals and the core states 
were recalculated after each iteration. We used 285, 240, 288 
and 455 inequivalent k points in the irreducible wedge of the 
body centered cubic (bcc), face centered cubic (fee), Llo and 
LI2 Brillouin zones, respectively. The random alloys (V-W, 
V-Fe) were treated by the coherent potential approximation. 19 
The equilibrium Wigner-Seitz radii (w) were extracted from 
the equation of state described by a Morse function 20 fitted to 
the total energies calculated for 17 different volumes around 
the equilibrium. All self-consistent calculations were carried 
out within the LDA and the gradient terms were included in 
the total energy within the perturbative approach. 21 The aver- 
age error (~ 7 x 10~ 4 Bohr in w) due to the above approxi- 
mation is below the numerical accuracy of our calculations. 22 

According to Table [I] neither PBE nor PBEsol yields sys- 
tematic errors for the Wigner-Seitz radius of transition met- 
als. Both of them underestimate the equilibrium volume of 3d 
metals (except for Cu) but overestimate those of the 5d met- 
als. As a consequence, the PBE/PBEsol over-binding remains 
nearly the same when going from pure Fe to Fe-V solid so- 
lutions and then to pure V (FigureQ] upper panel). However, 
for the V-W binary alloy, the inconsistent errors for V and W 
result in a significantly larger theoretical (PBE/PBEsol) vol- 
ume versus composition slope than the measured one (Figure 
[T] lower panel). 

Many of the recent density functionals have been adjusted 
to model systems^—i or to metal surfacesi— ^ Following the 
same idea, we investigate whether it is possible to adjust a 
GDF so that it produces systematic errors for the bulk prop- 
erties of metals. To this end, we consider the parameters /i 
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FIG. 1 . (Color online) Equilibrium Wigner-Seitz radii for Fe-V and 
V-W binary alloys. The present theoretical results (PBE/PBEsol: 
dashed lines, QNA: squares) are compare to Vegard's law (dashed- 
dotted line) and to experimental data (open and filled circles) i-*™ 



and f3 from the PBE/PBEsol scheme and compute the equi- 
librium volume of the selected mono -atomic metals as a func- 
tion of these "variables" for 0.075 < [i < 0.364 and 0.015 < 
f3 < 0.108. In figure [2] we show w(/i, (3) for V, Fe and Au 
relative to the corresponding experimental values (Table Q). 
We find that the errors in the theoretical w change smoothly 
from negative to positive values as a function of /i but are 
non-monotonous in terms of (3. As an consequence, changing 
Ppbe to /?pbeso1, for instance, results in minor effect on the 
theoretical volume of V or Fe. Another interesting feature is 
the positive (3 versus /i slope of the iso-error contour lines for 
P ^5 /^pbe/pbesoI- Hence, larger /i requires larger (3 to pre- 
serve the error in the equilibrium volume. This is the often 
quoted "error cancelation" between exchange and correlation 
terms. 

We find that for each element infinitely many pairs of 
{/U, /3} yield vanishing error in the equilibrium volume. These 
combinations form a continuous line in the {/x, /3}-space 
(marked by 0.000 in Figure The particular values /i opt 
corresponding to f3 opt = f3 PBE and (3 opt = /ipBEsoi are 
listed in Table [TT] We immediately observe that these "opti- 
mal" {/j,, /3} op t parameters are element specific. For instance, 
V and Fe require /i opt ~ 0.26 — 0.27, whereas for the 5d el- 
ements /i opt drops below ~ 0.13. These demonstrate that it 
is not possible to define a unique pair of "optimal" {/i, (3} pa- 
rameters. In other words, at least within the PBE/PBEsol con- 
straint, it is not possible to find a GDF that performs equally 
well for all metals. 

In an attempt to find the "best" {/x, /?} combination that 
leads to the smallest error in the equilibrium volume for the 
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FIG. 2. (Color online) Contours of relative errors in the Wigner-Seitz radius of V, Fe and Au. Specific sets of parameters corresponding to 
PBE, PBEsol, {/j,, P} opt with ft = /3pbe and /3 = /?pbeso1 are shown by symbols. The dotted grid represents those {/i, /?} pairs for which the 
calculations were performed. At every grid point, the equilibrium volume was derived from the total energy versus volume curves obtained 
using the particular {/i, /?} parameters. 



present set of mono-atomic solids, we minimized the cost 
function C(p,/3) = [u>,;(/i, /3)/iUi(expt) — l] 2 (i runs 
over the seven selected solids). Extending the search to < 
fx < 0.400 and < (3 < 0.260 resulted in ^ g _ opt = 0.151990 
and /3g_ op t = 0.230019 corresponding to C g _ op t = 0.0004. 
In terms of the cost function, the above "globally" opti- 
mized pair of parameters represents marginal improvement 
over C(PBE) = 0.0012 and C(PBEsol) = 0.0008. 

In order to reveal the origin of the difference be- 
tween the results obtained with PBE/PBEsol and those us- 
ing {n, /3} op t (leading to vanishing errors in w), we moni- 
tor the ratio between the corresponding enhancement func- 
tions F°^ l (r s , s)/Fj c BE (r s , s). Figure [3] displays the above 
ratio within the (001) plane of bcc V and fee Au. We no- 
tice that for both systems F^frsjs) and F^ E (r s ,s) are 
close to each other around the cell boundary but show large 
(positive for V and negative for Au) deviations for points in- 
side the atoms. That is, the interstitial region corresponds to 
nearly LDA regime, whereas the valence-core overlap region 



TABLE II. Special values for /x opt corresponding to /? op t = /3pbe 
and /3 op t = /3pbeso1 for the seven selected metals. 



Element 


/3pBE 


/3pBEsol 


V 


0.2722 


0.2646 


Fc 


0.2718 


0.2570 


Cu 


0.1732 


0.1587 


Nb 


0.1874 


0.1750 


Pd 


0.1284 


0.1124 


W 


0.1317 


0.1185 


Au 


0.1080 


0.0900 



is where the details of the GDF become important. This ob- 
servation is in line with those discussed by Fuchs et al.— and 
Csonka et al.— 

We conclude that for metals the primary error of a GDF 
is of local nature: the truly gradient-sensitive region is always 
localized around the atomic sites and the region in between the 
sites is less sensitive to the details of the density functional ap- 
proximation. This observation opens the door to an alternative 
SFA. According to that, for a solid the optimal subsystems are 
the individual (element specific) valence-core overlap regions 
which are connected by the nearly homogeneous valence elec- 
tron sea. For multi-component alloys, a possible realization of 
such quasi-non-local gradient-level approximation is the su- 
perposition of the component-optimized gradient-level func- 
tionals. Mathematically a QNA functional may be expressed 
as 

£ x Q c NA M = E / d 3 re™ A (n)FT« (r., (2) 

where Vl q represents the Wigner-Seitz cell around atom q 
and Fxc 9 is the PBE/PBEsol enhancement function based on 
{p,, /3} pt optimized for the alloy component q. Since for any 

{H,(3} 0?tq , F°c tq (r s ,s) reduces to - F^ A (r s ) within the 
interstitial region (characterized by s — > 0), to a good approx- 
imation the kernel of the functional from (|2]) is continuous at 
the cell boundaries. 

We demonstrate the above QNA scheme in the case of bcc 
V-Fe and V-W solid solutions and CuAu and CU3AU inter- 
metallic compounds adopting the Llo and LI2 structures, re- 
spectively. Using /3} op t obtained^ for bcc V, Fe and W, 
we computed the equilibrium Wigner-Seitz radii of V-Fe and 
V-W binary alloys as a function of chemical composition. The 
nearly perfect agreement between the present theoretical re- 
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FIG. 3. (Color online) Contour plot of the F°f (r s ,s) /F™ E (r s ,s) 
ratio within the (001) plane of the bcc V and fee Au. The atoms 
are located at the origin (x = y = 0), and the plots include points 
inside the sphere circumscribed to the Wigner-Seitz cells. The cross 
section between the cell boundary and the (001) plane is marked by 
thin solid lines. 



suits (denoted by QNA in Figure[T|) and the measured equilib- 
rium radii shows that the proposed component-optimized ap- 
proximation performs well for both dilute and concentrated al- 
loys. For CuAu and CU3AU the PBE/PBEsol errors in w rela- 
tive to the experimental values— are 1.6/0.3% and 1.1/0.7%, 
respectively. Using QNA with {/i, /3} op t optimized for fee Cu 
and Au, the above errors drop below 0.1%. 

It has been found that although for each elemental metal- 
lic solid there are several optimal GDFs that give vanish- 
ing errors in the theoretical equilibrium volume, no globally 
accurate GDF exists. Starting from the concept of subsys- 
tem functional approach, we have introduced a quasi-non- 
local gradient-level approximation based on the component- 
optimized GDFs. The so constructed QNA has been shown to 
perform equally well for mono-atomic and multi-component 
metallic systems. The proposed scheme can easily be imple- 
mented in any density functional method, including pseudo- 
potential methods via properly designed pseudo-potentials. 
Finally we should mention that QNA can be turned fully ab 
initio by using, instead the experimental equilibrium volumes, 
data provided by accurate many-body solvers. 
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